EXECUTIVE  SUMMARY 


Numerical  simulations  of  the  deformations  and  breakup  of  drops  have  been  done.  A  finite 
difference/front  tracking  numerical  technique  is  used  to  solve  the  unsteady  axisymmetric 
Navier-Stokes  equations  for  both  the  drops  and  the  ambient  fluid.  Inertial,  viscous,  and 
surface  tension  effects  are  all  accounted  for.  The  simulations  help  determine  where  in 
parameter  space  the  various  breakup  modes  take  place,  how  long  breakup  takes,  and  what 
the  resulting  drop  size  distribution  is.  The  goal  of  the  investigation  is  to  provide  results  that 
extend  and  complement  experimental  investigations,  and  lead  to  better  engineering  models 
of  drops  in  sprays. 

In  order  to  simulate  the  drop  breakup  at  small  density  ratios  (relevant  to  fuel  sprays  near 
critical  conditions),  two  cases  are  examined  in  detail:  a  density  ratio  close  to  unity 
( Prf/Po  =  1  •  15)  and  a  density  ratio  of  10.  In  both  cases,  the  drops  are  accelerated 
impulsively.  The  results  with  the  lower  density  ratio  can  be  re-scaled  to  predict  the  drop 
behaviors  at  other  density  ratios  close  to  one,  by  applying  the  Boussinesq  approximation. 

In  addition  to  full  simulations  where  the  Navier-Stokes  equations  are  solved,  a  few  inviscid 
simulations  have  also  been  done  to  identify  the  physical  mechanism  causing  drop  breakup. 


Four  non-dimensional  numbers  govern  the  breakup  of  drops.  In  addition  to  the  density  and 
the  viscosity  ratio,  the  ratio  of  inertia  to  surface  tension  is  described  by  the  Weber  number, 
whereas  the  ratio  of  inertia  to  viscous  force  is  represented  by  the  Reynolds  number. 
Sometimes,  these  two  numbers  are  combined  to  yield  the  Ohnesorge  number 
iPh^'^WelRe),  which  is  the  ratio  of  the  viscous  force  to  the  surface  tension.  The 
simulations  have  resulted  in  a  fairly  complete  picture  of  the  evolution  at  small  density 
ratios.  When  the  effect  of  viscosity  is  not  significant  {Oh  <  0.1),  the  Weber  number  has  the 
dominant  effects;  For  small  Weber  numbers,  the  drops  display  oscillatory  deformation, 
independent  of  the  density  ratio,  since  the  surface  tension  limits  the  deformation.  As  the 
Weber  number  is  increased,  the  drops  deform  into  a  disk-like  shape  first,  due  to  high 
pressure  at  the  fore  and  aft  stagnation  points  and  low  pressure  around  the  equator.  As  the 
deformation  progresses,  the  bulk  of  the  drop  fluid  ends  up  in  a  thick  rim  connected  by  a 
thin  film.  At  moderate  Weber  numbers,  the  rim  eventually  moves  faster  than  the  film  and  a 
backward-facing  bag  is  formed.  This  backward-facing  bag  mode  is  only  observed  for  the 
density  ratio  10  simulations.  In  the  lower  density  ratio,  1.15,  the  difference  between  the 
drop  and  the  fluid  inertia  vanishes  quickly  and  the  surface  tension  pulls  the  drops  back  into 
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L  BACKGROUND 


In  spray  combustion,  secondary  atomization  of  liquid  droplets  plays  an  important  role  in  the 
increase  of  surface  area  and  the  enhancement  of  heat  and  mass  transfer  between  the  fuel  and  the 
ambient  gas.  In  a  separate  report  provided  previously,  we  have  discussed  in  detail  the  secondary 
breakup  of  liquid  drops  that  are  accelerated  by  a  constant  body  force.  The  Eotos  number,  which  is 
the  ratio  of  the  body  force  to  the  surface  tension,  was  found  to  be  the  main  controlling  parameter. 
The  fluid  viscosities  and  the  density  ratio  had  secondary  effects.  Here,  we  focus  on  another  type 
of  disturbance  causing  drop  breakup:  impulsive  acceleration  of  the  droplet.  This  corresponds  to 
an  experimental  situation  where  a  liquid  droplet  is  accelerated  by  a  shock  wave.  This  type  of 
disturbance  was  employed  in  many  experimental  studies  of  secondary  breakup.  Earlier  studies 
include  Lane^,  Hinze^,  Hanson  et  al}.  Ranger  and  Nicholls^,  and  Gel’fand  et  al}.  For  more 
recent  studies  of  liquid  atomization,  see  Krzeczkowski®,  Borisov  et  alJ,  Reitz  and  Bracco®,  Pilch 
and  Erdman®,  Wierzba^®,  Hsiang  and  Faeth^^”^^,  and  Joseph  et  al}^. 

Experimental  results  are  usually  presented  in  terms  of  four  non-dimensional  numbers:  the  We¬ 
ber  number  We,  the  Reynolds  number  Re,  the  density  ratio  and  the  viscosity  ratio.  The  Reynolds 
number  is  sometimes  replaced  by  the  Ohnesorge  number.  These  numbers  are  defined  later  in 
the  paper.  The  experimental  data  shows  that  the  breakup  of  drops  in  gas  flows  can  generally  be 
categorized  into  six  different  modes^®: 

1.  "Vibrational  breakup  mode  where  the  original  drop  disintegrates  into  two  or  four  equal-sized 
smaller  drops. 

2.  Bag  breakup  mode  where  the  original  drop  deforms  into  a  torus-shaped  rim  spanned  by 
a  thin  fluid  film  that  ruptures  into  tiny  droplets,  followed  by  the  disintegration  of  the  rim 
into  larger  droplets.  Sometimes,  a  tiny  streamer  is  also  formed  inside  the  bag,  called  the 
parachute  or  the  umbrella  breakup  mode. 

3.  Transitional  mode  where  irregular  or  chaotic  fragmentation  takes  place.  In  this  mode,  the 
drop  breaks  up  due  to  both  the  formation  of  a  bag  and  tbe  stripping  of  fluid  from  the  drop 
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surface  in  the  form  of  stretched  fluid  filaments. 


4.  Shear  breakup  mode  where  continuous  stripping  of  fluid  from  the  drop  surface  occurs. 

5.  Explosive  breakup  mode  where  strong  surface  waves  disintegrate  the  drop  in  a  violent  man¬ 
ner. 

In  modes  2-4,  the  breakup  process  is  initiated  by  the  axisymmetric  flattening  of  the  drop  into  an 
ellipsoid  or  disk  shape.  For  modes  1  and  5,  the  formation  of  a  disk  is  not  seen  and  there  are  good 
reasons  to  believe  that  the  breakup  is  fully  three-dimensional  from  the  start.  This  categorization 
and  terminology  are  somewhat  arbitrary  and  different  variations  have  been  suggested  by  different 
researchers.  For  example,  mode  2  has  also  been  called  “bag-and-stamen  breakup”  by  Pilch  and 
Erdman®  and  “bag-jet  mechanism”  by  Krzeczkowski®. 

The  transitions  between  the  breakup  modes  described  above  take  place  at  critical  Weber  num¬ 
bers  that  generally  depend  only  weakly  on  the  Ohnesorge  number.  For  low  viscosity  drops,  the 
critical  Weber  numbers  are  approximately  10, 20, 30,  60,  and  1000,  respectively  for  the  transitions 
from  one  mode  to  the  next  one.  It  is  also  found  that  the  critical  Weber  numbers  increase  as  the 
Ohnesorge  number  increases.  Note  that  the  first  four  breakup  modes  are  realized  for  a  very  narrow 
range  of  Weber  numbers.  Wierzba^°  found  that  in  the  range  of  11  <  We  <  14,  the  drop  deforma¬ 
tion  and  the  breakup  are  very  sensitive  to  small  fluctuations  of  the  experimental  conditions  and  he 
could  further  refine  the  observed  breakup  into  five  sub-categories.  There  are  large  variations  in  the 
critical  Weber  numbers  in  the  available  experimental  data  due  to  different  test  conditions.  Thus 
the  numbers  presented  above  should  be  considered  only  as  a  rough  guide. 

Most  previous  studies  have  been  done  for  liquid-gas  systems  at  atmospheric  pressure  and  tem¬ 
perature.  The  behavior  of  liquid  droplets  at  smaller  density  ratios  (characteristic  of  high  pressure 
combustion  systems)  is  not  well-understood  yet.  In  this  study,  computational  results  are  presented 
for  impulsively  started  drops  in  which  the  motion  and  deformation  of  the  drop  is  caused  by  the 
initial  momentum  due  to  a  step  change  in  the  relative  velocity.  Numerical  simulations  have  been 
done  for  two  density  ratios,  1.15  and  10.  These  values  are  in  the  range  of  typical  density  ratios 
encountered  in  practical  spray  combustion  systems  (density  ratios  of  0(10)  for  diesel  engines^® 
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and  0(1)  for  rocket  motors  are  common).  Simulations  in  the  Boussinesq  limit  using  the  smaller 
density  ratio,  1.15,  can  be  extended  to  predict  the  drop  behavior  at  other  density  ratios  close  to 
1.  The  effects  of  the  individual  governing  parameters  are  examined  and  the  physical  mechanisms 
associated  with  different  breakup  modes  are  discussed. 


n.  FORMULATION  AND  NUMERICAL  METHOD 


A.  Navier-Stokes  equations 


The  physical  problem  and  the  computational  domain  are  sketched  in  Figure  1.  The  domain  is 
axisymmetric  and  the  left  boundary  is  the  axis  of  symmetry.  We  solve  for  the  motion  everywhere, 
both  inside  and  outside  the  drop.  The  Navier-Stokes  equations  are  valid  for  both  fluids,  and  a 
single  set  of  equations  can  be  written  for  the  whole  domain  as  long  as  the  jump  in  density  and 
viscosity  is  correctly  accounted  for  and  surface  tension  is  included. 

The  Navier-Stokes  equations  in  conservative  form  are: 

+  V  •  (puu)  =  -Vp  -I-  V  •  n{Vu  +  Vu’’)  +  cr  /  5(x  -  Xf)(n  x  V)  x  ndS.  (1) 
at  Js 


Here,  u  is  the  velocity,  p  is  the  pressure,  and  p  and  p,  are  the  discontinuous  density  and  viscosity 
fields,  respectively,  a  is  the  surface  tension  and  5  is  a  three-dimensional  delta  function.  Here,  the 
surface  tension  is  treated  as  a  singular  body  force  and  its  contribution  is  limited  to  the  interface 
itself.  The  integral  over  the  surface  of  the  drop,  S,  results  in  a  force  distribution  that  is  smooth  and 
continuous  along  the  drop  surface. 

The  above  equations  are  supplemented  by  the  incompressibility  condition: 


Vu  =  0 


(2) 


and  equations  of  state  for  the  physical  properties: 


2£  =  0;^  =  0 
Dt  ’  Dt 


(3) 


where  D/Dt  is  the  total  derivative.  These  two  equations  simply  state  that  the  physical  properties 
within  each  fluid  remain  constant  in  the  case  of  immiscible  fluids. 
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Dimensional  analysis  shows  that  a  set  of  four  independent  dimensionless  parameters  can  be 
constructed  for  the  dynamics  of  drop  deformation  and  breakup.  When  the  drop  is  subject  to  an 
impulsive  acceleration,  it  is  convenient  to  choose  the  Weber  number,  W e,  the  Reynolds  number 
based  on  the  ambient  fluid  properties.  Re,  density  ratio,  and  the  Reynolds  number  based  on  the 
drop  fluid  properties.  These  are  defined  by: 


...  PoU^D  ^  poUD  Pi  piVD 

We  =  - ;  Re  = - ;  —  ;  Rea  = - . 

po  P>d 


(4) 


Sometimes,  the  drop  Reynolds  number  is  replaced  by  viscosity  ratio,  PdfPo-  Here,  D  is  the  initial 
diameter  of  the  drop  and  U  is  the  initial  relative  velocity  between  the  drop  and  the  ambient  fluid. 
The  subscripts,  d  and  o,  denote  the  properties  of  the  drop  and  the  ambient  fluid,  respectively. 
Another  parameter  often  used  in  the  literature  is  the  Ohnesorge  number,  Oh  =  pf  ^JpDa  = 
y/We/Re,  which  is  the  ratio  of  the  viscous  shear  to  the  surface  tension. 

The  momentum  equations  and  the  continuity  equation  are  discretized  using  an  explicit  second- 
order  predictor-corrector  time-integration  method  combined  with  a  second-order  centered  differ¬ 
ence  method  for  the  spatial  derivatives.  The  discretized  equations  are  solved  on  a  fixed,  staggered 
grid,  order  to  improve  the  accuracy  of  the  computation  in  some  specific  regions,  non-uniform 
grid  spacing  is  used.  The  full-slip  boundary  condition  is  applied  to  all  four  boundaries  unless  spec¬ 
ified  otherwise.  For  simulations  with  pd/ Po  =  1-15,  a  fixed  computational  domain  is  used  whereas 
for  the  pd/po  =  10  simulations,  a  reference  frame  moving  with  the  centroid  of  the  drop  is  used. 
This  enables  us  to  keep  the  computation  at  a  reasonable  cost,  while  maintaining  the  necessary 
resolution  around  the  moving  drop.  In  order  to  accurately  simulate  a  highly-deformed  interface 
during  breakup,  a  front-tracking  method  is  employed,  where  the  interface  is  represented  by  marker 
points  that  are  moved  by  interpolating  their  velocities  from  the  stationary  grid.  These  points  are 
connected  to  form  a  “front”  which  is  used  to  construct  the  property  fields  at  each  time  step.  This 
front  is  also  used  to  calculate  the  surface  tension.  For  a  more  detailed  description  of  the  front 
tracking  method,  see  Tryggvason,  Bunner,  Ebrat,  and  Tauber^®.  The  majority  of  the  simulations 
presented  here  were  carried  out  on  HP  9000  workstations.  A  typical  run  generally  required  be¬ 
tween  4000  and  120000  time-steps  and  took  12-240  hours,  depending  on  the  parameters  of  the 
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problem. 


B.  Inviscid  flow 

To  address  to  what  extend  the  drop  evolution  can  be  described  by  an  inviscid  model,  a  few 
simulations  were  done  using  a  vortex  method.  The  interface  separating  the  drop  and  the  ambient 
fluid  is  a  vortex  sheet  and  an  evolution  equation  for  the  vortex  sheet  strength  7  =  (u^  -  u<,)  •  s 
can  be  derived  by  subtracting  the  tangential  components  of  the  Euler  equations  on  either  side  of 
the  interface.  The  resulting  equation  is  (see,  for  example,  Tryggvason,  1988): 
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Dt  ^  ds 
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\di^ 

s  +  --7^ 


1  dj^  \  a  Ok 

%  ds'  Pd  +  Pods' 


Here,  A  =  {pd  -  po)l{pd  +  Po)  is  the  Atwood  number,  U  =  (u^  +  Uo)/2  is  the  average  of 
the  velocities  on  either  side  of  the  vortex  sheet,  and  k  is  the  mean  curvature.  Given  the  vortex 
sheet  strength  7,  the  velocity  is  foimd  by  the  Biot-Savart  integral.  For  computational  purpose, 
the  axisymmetric  vortex  sheet  is  discretized  by  a  finite  number  of  vortex  rings.  The  azimuthal 
integration  can  be  done  analytically  and  the  integral  is  therefore  replaced  by  a  summation  over  the 
discrete  vortex  rings.  The  radial  and  axial  velocities  at  a  point  on  ring  j  are  given  by 

%  =  ^  I  [(fe|)  m>) + (6) 


Tinkij  ^{{ri  +  rj)kl-2ri 


”  TT  ^  {4rj 


1  - 


E{kii)A2riK{kii)\ 


Here  K{k)  is  the  Complete  Elliptic  integral  of  the  first  kind,  and  E{k)  is  the  Complete  Elliptic 


integral  of  the  second  kind  and 
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The  Elliptic  integrals  can  be  computed  efficiently  by  a  polynomial  approximation^^.  When  the 
axisynunetric  vortex  sheet  is  replaced  by  discrete  vortex  rings,  the  rings  must  be  given  a  finite  core 
size  to  avoid  infinite  self-induced  velocity.  This  can  be  accomplished  simply  by  replacing  kij  by 
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where  5  is  a  small  regularization  parameter  that  determines  the  size  of  a  “vortex  blob”.  In  the  limit 
of  iV  00  and  5  ^  0  the  solution  will  be  independent  of  the  exact  value  of  5  (except  at  isolated 
points  where  roll-up  takes  place). 

The  evolution  equation  for  the  circulation  of  each  vortex  ring  can  be  obtained  by  integrating 
Equation  (5)  over  a  small  material  segment: 
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where  a"''  and  a~  are  the  two  end  points  of  the  segment.  Equation  (10),  along  with  the  equation 
D:x./Dt  =  U,  where  x  =  (r,  z)  is  the  position  vector  for  each  vortex  ring,  can  be  integrated  using 
a  time-integration  method  such  as  the  fourth-order  Runge-Kutta  method,  once  the  initial  vorticity 
distribution  and  the  geometry  are  given.  The  DV/Dt  term  is  found  in  an  iterative  manner:  First 
the  interface  is  advanced  and  the  circulations  are  updated  using  DTJ/Dt  from  last  time-step  as  an 
approximation.  Then  the  velocities  at  the  new  position  are  found  using  the  approximate  value  of 
the  circulations.  The  velocities  are  then  used  to  calculate  DU / Dt  at  the  original  interface  position 


and  the  process  is  repeated. 


m.  SECONDARY  BREAKUP  OF  DROPS  DUE  TO  IMPULSIVE  ACCELERATION 

Most  of  the  simulations  presented  here  are  for  pd/ Po  =  10.  The  effects  of  the  Weber  number, 
the  Reynolds  number,  and  the  viscosity  ratio  are  studied.  In  addition,  a  few  simulations  have 
been  done  for  Pd/Po  =  1-15  to  examine  the  effect  of  the  density  ratio.  As  discussed  in  part  I, 
the  Boussinesq  approximation  applies  when  the  drop  and  the  ambient  fluid  have  similar  densities. 
Thus  a  single  simulation  with  one  specific  value  of  the  density  ratio  close  to  1  can  be  rescaled 
for  a  range  of  density  ratios.  When  presenting  the  results,  time  is  non-dimensionalized  with  the 
diameter  and  the  initial  relative  velocity,  t*  =  t/{D/U). 
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A.  Effect  of  We  at  a  fixed  Re 


1.  Evolution  of  drops  with  Pd/Po  —  10 

In  Figure  2,  the  deformation  and  breakup  of  drops  with  pd/po  =  10  are  shown  for  Re  =  242 
and  Red  =  1935.  Results  for  We  =  3.74,  12.5,  18.7,  28.1,  37.4,  46.8,  56.1,  74.8,  93.5  and  oo 
are  compared.  For  these  parameters  (excluding  the  We  =  oo  case),  the  Ohnesorge  numbers.  Oh, 
are  between  0.008-0.04.  Experimental  evidence  suggest  that  when  Oh  <  0.1,  the  viscosity  has 
only  minor  effects  and  that  surface  tension  effects  are  dominant.  The  computations  were  done 
using  a  resolution  of  256  x  512  grid  points.  The  drop  in  (a)  with  We  =  3.74  shows  oscillatory 
deformation.  The  drops  with  We  =  12.5  and  18.7  in  (b)  and  (c)  initially  develop  an  indentation  at 
the  top  and  it  continues  to  deepen.  As  the  deformations  increase  the  drag,  the  initial  momentum  of 
the  drop  fluid  decreases.  This  reduces  the  aerodynamic  pressure  variations  around  the  drop  and  the 
influence  of  the  surface  tension  starts  to  increase.  At  later  stages  of  the  deformation,  a  backward- 
facing  bag  is  formed.  Eventually,  the  backward-facing  bags  in  (b)  and  (c)  stop  growing  and  the 
drops  restore  their  originally  spherical  shapes.  The  drop  shown  in  (d)  for  We  =  28.1  shows  the 
formation  of  a  backward-facing  bag  more  clearly.  The  Weber  numbers  for  the  drops  in  (b)-(d) 
are  in  the  range  where  a  bag  breakup  mode  is  observed  experimentally  for  higher  density  ratios. 
The  results  in  (e)-(i),  for  even  higher  Weber  numbers,  show,  on  the  other  hand,  the  formation  of  a 
forward-facing  bag.  The  drop  with  zero  surface  tension  in  (j)  also  displays  a  forward-facing  bag. 
In  this  case  small  scale  irregularities  are  observed  both  on  the  edge  and  on  the  top  surface  of  the 
drop.  This  is  due  to  the  absence  of  the  surface  tension,  which  suppresses  the  formation  of  the 
small-scale  structures  caused  by  the  interfacial  shear. 

The  aspect  ratio  (the  maximum  width  of  the  drop  divided  by  its  thickness  at  the  centerline), 
centroid  velocity,  and  surface  area  of  some  of  the  drops  shown  in  Figure  2  are  plotted  versus  non- 
dimensional  time,  t*,  in  Figure  3.  The  aspect  ratio  plot  displays  an  oscillation  of  the  drop  with  the 
lowest  We  and  monotonic  decrease  for  the  other  drops.  When  the  drops  with  We  =  12.5  and  18.7 
start  to  resume  their  initial  shape,  an  abrupt  increase  in  the  aspect  ratio  is  observed.  The  velocity 
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plot  shows  a  monotonic  decline  for  all  W^e  shown.  The  surface  area  plot  shows  that  the  rate  of 
deformation  increases  with  We. 

In  Figure  4  and  5,  vorticity  contours  (left)  and  streamlines  (right)  at  selected  times  are  plotted 
along  with  the  drop  contour  for  the  We  =  28.1  and  We  =  93.5  drops  shown  in  Figure  2(d)  and 
(i)  respectively.  In  both  cases,  the  vorticity  plots  show  that  most  of  the  vorticity  is  created  at  the 
outer  edge  of  the  drop,  as  expected.  The  streamline  pattern  for  the  We  =  28.1  drop  shows  that  the 
backward-facing  bag  is  stretched  upward  when  the  wake  downstream  of  the  drop  detaches  from 
it.  On  the  other  hand,  the  closed  streamlines  around  the  We  =  93.5  drop  suggests  that  the  drop 
moves  as  a  vortex  ring,  forming  a  forward-facing  bag. 

Inviscid  flow  simulations  using  a  vortex  method  are  presented  for  We  =  93.5  and  pd/ po  =  10 
in  Figure  6.  The  evolution  of  three  inviscid  drops  with  different  blob  sizes,  5  =  0.0125, 0.01,  and 
0.0075,  are  shown  in  (a),  (b),  and  (c),  respectively.  In  addition,  the  viscous  simulation  for  the  same 
Weber  number  and  the  density  ratio  — ^already  shown  in  Figure  2(i)) —  is  illustrated  again  in  (d) 
for  comparison.  In  each  column,  the  drop  shapes  are  plotted  every  t*  —  0.4835.  Note  that  5  was 
introduced  in  Equation  (9)  only  to  avoid  singularity  and  that  the  choice  of  5  is  somewhat  arbitrary. 
The  overall  drop  evolution  is,  however,  not  very  sensitive  to  the  exact  value  of  5:  In  all  three  cases, 
a  thin  film  of  drop  liquid  is  pulled  away  from  the  shoulder  of  the  main  drop.  Differences  exist  in 
the  small  scale  structures  such  as  the  thickness  of  the  film,  which  reduces  faster  for  smaller  6.  The 
viscous  drop  in  (d)  also  displays  an  evolution  similar  to  those  of  the  inviscid  drops. 

The  centroid  velocities  of  the  drops  in  Figure  6  (a)-(d)  are  compared  in  Figure  7.  As  the  drops 
starts  to  move,  the  viscous  simulation  shows  a  faster  decrease  of  the  initial  velocity,  compared  to 
the  inviscid  results.  Later,  the  inviscid  drops  deform  more  and  the  trend  is  reversed.  The  degree  of 
deformation  of  the  inviscid  drops  also  increases  weakly  as  S  becomes  smaller,  so  the  velocity  of 
an  inviscid  drop  with  a  smaller  6  decreases  faster. 

In  spite  of  some  differences,  due  to  viscosity  on  the  one  hand  and  finite  regularization  on  the 
other,  the  viscous  drop  and  the  inviscid  ones  display  essentially  the  same  behavior.  Traditional 
“shear”^  or  “boundary-layer  stripping”'^  breakup  theories  are  based  on  the  assumption  that  viscous 
boundary  layers  develop  on  both  side  of  the  interface  and  viscous  shear  causes  a  portion  of  the  drop 
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fluid  to  be  stripped  from  the  drop  shoulder.  The  validity  of  this  assumption  has  been  questioned  by 
Lee  and  Reitz^^,  whose  experimental  observations  in  the  shear  breakup  regime  were  independent 
of  the  Reynolds  number,  suggesting  weak  influence  from  the  viscous  boundary  layers.  As  an 
alternative,  they  proposed  a  “stretching-thinning”  type  breakup  mechanism  where  the  thin  edge 
of  the  drop,  flattened  by  the  aerodynamic  pressure,  is  deflected  in  the  downstream  direction  by 
blowing  of  the  ambient  fluid.  Our  inviscid  simulations  show  that  the  initial  formation  of  the  thin 
film  from  the  drop  shoulder  is  an  inviscid  phenomenon,  supporting  the  “stretching-thinning  type” 
mechanism. 


2.  Evolution  ofdmps  with  Pd/Po  =  1*15 

Figure  8  shows  the  deformation  versus  time  of  drops  with  a  density  ratio  of  1.15  in  the  same 
way  as  in  Figure  2.  The  Weber  numbers  are  2.73,  13.7,  27.4,  54.7,  and  oo,  the  Reynolds  number 
is  fixed  at  331,  and  the  drop-based  Reynolds  number  is  381 .  The  computations  were  done  using  a 
128  X  256  grid.  The  low  Weber  number  drop  in  (a)  oscillates  due  to  the  high  surface  tension.  As 
We  is  increased  to  13.7,  the  drop  starts  to  develop  an  indentation  at  the  top  as  shown  in  (b)  but 
as  it  falls,  the  momentum  of  the  drop  decreases  and  the  surface  tension  causes  it  to  oscillate.  The 
drop  shown  in  (c)  for  We  =  27.4  also  deforms  into  an  indented  ellipsoid  initially.  The  indentation 
deepens  progressively  and  later  meets  the  bottom  of  the  drop  surface,  forming  a  forward-facing 
bag  as  observed  in  the  continuous  acceleration  case.  However,  since  there  is  no  force  to  maintain 
the  motion,  the  drop  eventually  resumes  its  initial  shape.  The  drop  with  We  —  54.7  in  (d),  on  the 
other  hand,  shows  increased  initial  deformation  which  results  in  a  bigger  rim  that  is  connected  by 
a  forward-facing  bag.  For  zero  surface  tension,  (e),  the  drop  displays  a  roll-up  of  the  interface. 
A  similar  roll-up  has  been  observed  in  the  constant  acceleration  cases  for  drops  with  no  surface 
tension. 

In  Figure  9,  the  normalized  aspect  ratio,  centroid  velocity,  and  surface  area  are  plotted  versus 
t*  in  (a)-(c),  respectively,  for  the  drops  shown  in  Figure  8.  The  aspect  ratio  plot  shows  shape  os¬ 
cillation  for  the  two  lower  values  of  W e.  For  higher  W e,  the  aspect  ratio  decreases  monotonically 
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to  zero  as  the  indentation  deepens  progressively.  The  velocity  plot  in  (b)  shows  a  rapid  decline 
initially.  Later,  the  velocity  of  the  drops  with  We  =  2.73  and  13.7  decreases  but  with  fluctua¬ 
tions.  The  velocities  of  the  drops  with  three  higher  Weber  numbers  decrease  monotonically  until 
it  reaches  a  minimum  and  starts  to  increase  again.  The  surface  area  plot  in  (c)  show  the  effect  of 
We  on  deformation:  As  We  increases,  the  slope  of  the  curve  increases. 

B.  Effect  of  Re 

In  order  to  see  the  effect  of  the  viscosity  of  the  ambient  fluid,  simulations  with  four  different 
Reynolds  numbers.  Re  =  387,  242,  121,  60.5  (from  left  to  right)  are  compared  in  Figure  10.  In 
the  top  row  (a).  We  =  18.7,  Pd/ Po  —  Red  —  1935.  The  overall  evolution  of  the  drops 

with  two  higher  Reynolds  numbers.  Re  =  387  and  242,  are  very  similar  and  the  only  difference  is 
in  the  small  scale  structure  of  the  drop  contour.  The  results  for  the  two  lower  Reynolds  numbers. 
Re  =  121  and  60.5,  however,  display  some  differences:  Increasing  the  ambient  fluid  viscosity 
(decreasing  the  Reynolds  number)  while  the  other  parameters  are  fixed  can  lead  to  a  transition 
from  the  backward-facing  bag  mode  to  the  oscillatory  deformation  mode.  Another  comparison  is 
made  in  (b)  for  We  =  74.8,  while  other  parameters  remain  the  same  as  in  (a).  Here,  a  transition 
from  the  forward-facing  bag  mode  to  the  backward-facing  bag  mode  is  observed  as  the  Reynolds 
number  is  reduced  from  121  to  60.5. 

Based  on  these  observations,  it  is  clear  that  progressively  higher  Weber  numbers  are  necessary 
in  order  to  observe  the  same  mode  of  deformation,  as  the  Reynolds  number  is  decreased.  The 
translation  of  the  boundaries  between  different  deformation  modes — oscillation,  backward-facing 
bag,  and  forward-facing  bag — ^to  higher  Weis  clearly  due  to  the  increased  viscous  dissipation. 

In  Figure  11,  the  surface  areas  of  the  drops  shown  in  Figure  10  are  plotted  versus  the  non- 
dimensional  time,  f*.  The  comparison  in  (a)  for  drops  with  We  =  18.7  indicates  that  the  rate  of 
increase  of  the  surface  area  becomes  higher  as  the  Reynolds  number  increases,  due  to  less  viscous 
dissipation.  The  same  trend  is  confirmed  for  We  =  74.8  in  (b). 
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C.  Effect  of  Red 


In  Figure  12,  the  effect  of  drop  viscosity  is  shown  for  drops  with  Pd/Po  —  10  and  Re  =  242. 
The  figures  in  each  frame  denotes  the  dimensionless  time  when  the  drop  is  plotted.  In  the  top  row, 
four  cases  are  compared  for  different  drop  viscosities  (represented  by  the  Reynolds  number  based 
on  drop  properties)  at  a  fixed  Weber  number.  We  =  28.1.  As  the  drop  viscosity  is  increased  from 
left  to  right  by  an  order  of  10^,  the  drop  deformation  is  greatly  reduced.  The  least  viscous  drop 
with  Red  =  1935  clearly  shows  the  formation  of  a  backward-facing  bag,  while  the  most  viscous 
drop  with  Red  =  1.935  remains  in  an  oblate  shape.  In  the  middle  row,  a  similar  comparison  is 
made  fot  We  =  56.1.  Again,  by  increasing  the  drop  viscosity,  the  drop  deformation  is  reduced 
and  the  drop  changes  from  a  forward-facing  bag  to  an  oblate  drop.  The  result  for  We  =  93.5 
drops  shown  in  the  bottom  row  also  displays  a  similar  trend. 

D.  Deformation  and  breakup  regime  map 

In  Figure  13(a),  a  breakup  map  is  shown  to  summarize  the  deformation  and  breakup  modes  of 
drops  with  Pd/Po  =  10  and  Red  =  1935.  The  horizontal  and  vertical  axes  represent  the  Reynolds 
number  and  the  Weber  number  based  on  the  ambient  fluid,  respectively.  The  various  breakup 
modes  are  denoted  by  different  symbols.  When  the  Weber  number  is  low,  surface  tension  pre¬ 
vents  large  deformation  so  the  drop  only  oscillates.  As  the  Weber  number  increases  past  a  critical 
value  (approximately  16)  for  high  Re,  a  backward-facing  bag  starts  to  emerge.  When  the  Weber 
number  is  higher  than  approximately  30,  the  drops  break  up  in  a  backward-facing  bag  mode  for 
low  Reynolds  numbers  and  in  a  forward-facing  bag  mode  at  high  Reynolds  numbers.  The  transi¬ 
tions  from  the  backward-facing  bag  mode  to  the  forward-facing  bag  mode  occur  at  progressively 
lower  Weber  numbers  as  the  Reynolds  number  increases.  As  the  Weber  number  is  raised  to  100, 
approximately,  all  drops  break  up  in  the  forward-facing  bag  mode  for  the  Reynolds  number  range 
examined.  Finally,  when  the  surface  tension  vanishes  (We  oo),  the  strong  shear  due  to  the 
outside  flow  peels  off  the  drop  interface  and  shear  breakup  is  observed  when  the  Reynolds  number 
is  greater  than  100. 


A  similar  breakup  map  is  presented  in  Figure  13(b)  for  the  drops  with  =  and  Rea,  = 
1935.  Again,  when  the  Weber  number  is  low,  the  drops  only  oscillate  due  to  large  surface  tension. 
As  the  Weber  number  is  increased,  the  backward-facing  bag  mode  displayed  in  (a)  for  paf  Po  =  10 
is  no  longer  observed  as  the  momentum  difference  between  the  drop  and  the  ambient  fluid  is 
relatively  smaller,  in  case  of  density  ratio  close  to  one.  Instead,  the  drops  in  this  Weber  number 
range  oscillate  back  into  their  initial  spherical  shapes.  When  the  Weber  number  is  greater  than 
100,  the  drop  deforms  into  a  forward-facing  bag  and  the  roll-up  of  the  interface  is  observed.  The 
effect  of  the  Reynolds  number  is  similar  to  what  has  been  observed  in  (a):  the  boundaries  between 
different  deformation/breakup  modes  move  to  higher  Weber  numbers  as  the  Reynolds  number 
decreases. 


IV.  SIMULATIONS  OF  HEAT  TRANSFER  DURING  BREAKUP 

The  ultimate  reason  for  atomizing  a  fuel  jet  is  to  increase  the  evaporation  rate  of  the  fuel.  Be¬ 
fore  evaporation,  the  fuel  must  heat  up.  As  fuel  drops  break  up,  their  surface  area  increases  greatly 
and  therefore  the  heat  transfer  to  the  drop.  Simulating  heat  and  mass  transfer  at  intermediate 
Reynolds  numbers,  in  addition  to  fluid  flows  is  a  very  complicated  problem.  Under  these  circum¬ 
stances,  various  transport  mechanisms  are  involved  and  the  accurate  determination  of  the  relative 
importance  of  each  mechanism  becomes  a  challenging  task.  Some  simplifications  are  possible 
when  the  drop  has  transport  properties  considerably  different  from  those  of  the  surrounding  fluid. 
For  problems  not  involving  a  phase  change,  a  comprehensive  discussion  of  these  approximations 
is  give  by  Abramzon  and  Borde^^.  When  both  fluids  have  comparable  transport  properties,  the 
problem  becomes  more  complex  and  fiill  numerical  simulations  are  often  employed.  For  example, 
studies  of  isolated,  as  well  as  interacting,  vaporizing  fuel  droplets  at  intermediate  Reynolds  num¬ 
bers  are  described  by  Chiang^®,  Chiang  and  Sirignano^^  and  Haywood  et  al.'^.  In  these  studies, 
simulations  are  limited  to  cases  with  either  small  or  moderate  deformation. 

Here,  we  examined  the  heat  transfer  during  secondary  breakup  of  drops  with  a  small  density 
ratio,  PdfPo  =  1-15.  The  drops  are  accelerated  by  a  body  force  and  have  an  initial  temperature 
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lower  than  the  surrounding  fluid.  Phase  change  and  density  variation  with  temperature  are  not 
considered.  The  temperature  field  at  each  time  step  is  obtained  by  solving  the  unsteady  energy 
equation.  If  heat  generation  due  to  viscous  friction  and  radiation  heat  transfer  is  neglected,  the 
energy  equation  is: 

+  V .  {pc^T)  =  V  •  kVT.  (1 1) 

at 

This  equation  is  solved  on  a  fixed  grid  by  an  explicit  second  order  method  in  the  same  way  as  the 
momentum  equation. 

In  Figure  14,  the  evolution  of  the  drop  interface  is  shown  for  different  Eo.  In  all  cases, 
p^/po  —  1.15.  Values  of  Ohnesorge  numbers.  Oho  =  0  05  and  Oha  =  0.0466,  are  chosen  so 
that  viscous  effects  are  small  compared  to  surface  tension  effects.  The  computations  were  done 
using  a  fixed  coordinate  system.  When  Eo  is  small,  the  drop  is  deformed  into  an  oblate  ellip¬ 
soidal  shape  and  reaches  a  steady  state  as  shown  in  (a)  for  Eo  =  2.4.  The  drop  shown  in  (b) 
for  Eo  =  24  initially  develops  an  indentation  at  the  top  and  later  deforms  into  a  disk-like  shape. 
Then  the  thickness  of  the  drop  near  the  symmetry  axis  continues  to  decrease  and  most  of  the  drop 
fluid  moves  outward  toward  the  edge  of  the  drop.  Finally,  the  center  portion  of  the  front  surface  is 
pushed  upward,  forming  a  backward-facing  bag.  At  this  stage,  most  of  the  drop  fluid  is  contained 
in  the  annular-shaped  rim.  As  time  progresses,  the  bag  portion  of  the  drop,  which  forms  a  thin 
membrane  spanning  the  rim,  expands  both  radially  outward  and  vertically  upward.  Experimental 
evidence  indicate  that  it  will  eventually  break  into  smaller  drops.  In  (c),  where  Eo  is  48,  the  drop 
remains  in  a  convex  shape  showing  no  sign  of  forming  a  bag  while  the  shape  oscillation  continues 
during  the  entire  time  simulated.  When  Eo  is  144  in  (d),  the  indentation  developed  at  the  top  of 
the  drop  deepens  continuously  until  it  reaches  the  bottom  of  the  drop,  creating  a  forward-facing 
bag.  Eventually,  the  heavier  edge  falls  faster  than  the  thin  bag. 

Figure  15(a)  shows  vorticity  contours  (left)  and  streamlines  (right)  when  the  drop  shown  in 
Figure  14(b)  develops  a  backward-facing  bag.  In  the  vorticity  plot,  most  of  the  vorticity  is  created 
at  the  outer  edge  of  the  drop,  as  expected.  The  streamlines  plotted  with  respect  to  a  reference  frame 
moving  with  the  drop  show  that  the  backward-facing  bag  is  stretched  upward  in  the  downstream 
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direction.  The  existence  of  the  wake  behind  the  drop  is  evident  when  the  backward-facing  bag 
is  formed.  Figures  15(b)  illustrates  vorticity  contours  (left)  and  streamlines  (right)  for  the  drop 
shown  in  Figure  14(d)  when  a  forward-facing  bag  is  formed.  Again,  much  of  the  vorticity  is 
created  at  the  outer  edge  of  the  drop.  However,  the  closed  streamlines  around  the  drop  suggest  that 
the  drop  is  essentially  a  vortex  ring.  The  size  of  the  wake  is  smaller  than  in  Figure  15(a). 

In  Figures  16,  the  surface  area  of  the  drop  S  (normalized  by  the  initial  value  So)  is  plotted 
versus  t*  for  the  drops  shown  in  Figures  14  (a)-(d).  Note  that  the  surface  area  increases  rapidly 
upon  drop  breakup.  For  example,  the  Eo  =  24  drop,  which  deforms  into  a  backward-facing  bag, 
shows  a  nearly  steady  increase  with  some  fluctuations.  The  slope  of  the  graph  is  steeper  when 
the  drop  is  about  to  form  the  backward-facing  bag  but  once  the  bag  is  formed  and  starts  to  grow 
in  size,  the  slope  is  reduced.  Meanwhile,  the  result  for  Eo  =  48  shows  an  oscillatory  behavior. 
The  magnitude  of  the  surface  area  fluctuations,  however,  does  not  grow.  The  graph  for  the  drop 
with  higher  144  illustrates  a  rapid  initial  increase  as  the  drops  deform  into  a  forward-facing  bag. 
It  then  decreases  when  the  draining  of  the  drop  fluid  from  the  symmetry  axis  is  completed.  When 
the  heavy  outer  rim  of  the  drop  containing  most  of  the  drop  fluid  falls  faster  than  the  thin  film 
near  the  symmetry  axis,  the  total  surface  area  of  the  drop  starts  to  rapidly  increase  again  due  to  the 
formation  of  a  backward-facing  bag. 

The  temperature  field  and  the  drop  contour  are  plotted  in  Figure  17  for  the  drops  shown  in 
Figure  14  (a),  (b),  and  (d).  In  all  cases,  the  drops  have  identical  thermal  properties,  Pr  =  1, 
kd/ko  =  0.2  and  Cp^d/cp,o  =  1-  The  spherical  drop  is  plotted  at  f  =  28,  the  bag-breakup  drop  at 
t  =  28,  and  the  shear  breakup  drop  at  f  =  4.  The  temperature  field  around  the  spherical  drop 
shows  that  recirculation  inside  the  drop  leads  to  a  ring  of  cold  fluid  inside  the  drop  and  a  steep 
temperature  gradient  near  the  front  of  the  drop.  As  the  drop  heats  up,  the  ambient  fluid  cools  down 
and  a  cold  thermal  wake  extends  downstream  from  the  drop.  The  temperature  field  inside  the  drop 
undergoing  bag  breakup  is  very  different.  There  is  essentially  no  recirculation  inside  the  drop  and 
the  coldest  fluid  is  at  the  centerline  of  the  drop.  Since  there  is  a  large  wake  behind  the  drop,  the 
cold  thermal  wake  is  much  larger  than  for  the  low  Eo  drop.  The  drop  undergoing  shear  breakup 
has  only  a  small  thermal  wake,  but  the  large  deformation  of  the  drop  leads  to  a  large  heat  transfer 
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rate. 


In  Figure  1 8,  the  evolution  of  the  average  temperature  is  presented  for  the  three  drops  examined 
above.  Initially  the  drops  gain  heat  at  approximately  the  same  rate,  but  as  they  start  to  deform, 
the  rate  depends  strongly  on  the  breakup  mode.  The  heat  transfer  rate  is  slowest  for  the  nearly 
spherical  drop  and  quickly  increases  as  the  Eo  increases.  As  Figure  16  showed,  the  rate  of  increase 
of  surface  area  increases  greatly  as  Eo  is  increased.  However,  the  flow  field  is  also  different  and 
convective  heat  transfer  changes.  How  much  of  the  increased  heat  transfer  can  be  attributed  to 
increased  surface  area  and  how  much  to  increased  convection  has  not  been  examined  in  detail  yet. 
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FIGURES 


FIG.  1.  Schematic  diagram  of  the  computational  domain. 

FIG.  2.  Evolution  of  impulsively  started  drops  with  Re  =  242,  Red  =  1935,  and  PdfPo  —  10.  Results 
for  ten  W^e  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  next  to  the  drops  denote 
t*  when  the  interfaces  are  plotted.  The  centroid  of  the  drops  in  a  column  are  separated  by  a  fixed  distance. 
The  gap  between  two  successive  drops  in  a  column  does  not  represents  the  distance  the  drop  travels  during 
the  time  interval.  The  computations  were  done  on  a  256  x  512  grid. 

FIG.  3.  Non-dimensionalized  aspect  ratio,  centroid  velocity,  and  surface  area  plotted  versus  <*  for 
selected  cases  of  the  drops  shown  in  Figure  2.  Results  for  five  PFe  =  3.74,  12.5, 18.7,  37.4,  and  93.5  are 
compared.  Re  =  242,  Red  =  1935,  and  pdfPo  =  10. 

FIG.  4.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  shown  in  Figure  2(d). 

FIG.  5.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drop  shown  in  Figure  2(i). 

FIG.  6.  Inviscid  flow  simulations  for  PFe  =  93.5  and  Pd/Po  =  10  with  three  different  values  of 
5  =  0.0125  (a);  0.01  (b);  0.0075  (c).  In  (d),  the  viscous  result  shown  in  Figure  2(i)  is  included  for  compari¬ 
son.  The  drop  shapes  are  plotted  every  t*  =  0.4835  in  each  column. 

FIG.  7.  Centroid  velocities  of  the  drops  shown  in  Figure  6  plotted  versus  t*. 

FIG.  8.  Evolution  of  impulsively  started  drops  with  Re  =  331,  Red  =  381,  and  Pd/Po  =  1-15.  Results 
for  five  PFe  are  shown  as  denoted  by  the  numbers  below  the  figures.  The  numbers  inside  the  frames  denote 
t*  when  the  interfaces  are  plotted.  The  computations  were  done  on  a  128  x  256  grid. 

FIG.  9.  Non-dimensionalized  aspect  ratio,  centroid  velocity,  and  surface  area  plotted  versus  i*  for  the 
drops  shown  in  Figure  8.  Results  for  five  PFe  =  2.73, 13.7,  27.4,  54.7,  and  oo  are  compared.  Re  =  331, 
Red  =  1935,  and  pd/Po  —  1-15. 
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FIG.  10.  Effects  of  the  Reynolds  number  are  shown  for  drops  with  pd/po  =  10  and  Red  =  1935.  In 
the  top  row  (a),  We  =  18.7  and  evolutions  of  drops  with  four  different  Reynolds  numbers,  Re  =  387, 
242, 121,  60.5  (from  left  to  right)  are  compared.  In  the  bottom  row  (b),  a  similar  comparison  is  made  for 
We  =  74.8  with  the  same  set  of  the  Reynolds  numbers.  The  numbers  next  to  the  drops  denote  t*  when 
the  interfaces  are  plotted.  The  centroid  of  the  drops  in  a  column  are  separated  by  a  fixed  distance.  The  gap 
between  two  successive  drops  in  a  column  does  not  represents  the  distance  the  drop  travels  during  the  time 
interval.  The  computations  were  done  on  a  256  x  512  grid. 

FIG.  11.  Non-dimensional  surface  area  plotted  versus  <*  are  compared  for  four  Reynolds  numbers. 
Re  =  60.5, 121, 242,  and  387.  In  (a),  results  for  the  We  =  18.7  drops  shown  in  Figure  10(a)  are  presented 
and  in  (b),  those  for  the  We  =  74.8  drops  in  Figure  10(b)  are  shown. 

FIG.  12.  Drop  viscosity  effect  on  the  deformation  of  drops  with  Pd/Po  =  10  and  Re  =  242.  The 
number  in  each  frame  denotes  the  dimensionless  time  when  the  drop  is  plotted. 

FIG.  13.  Breakup  mode  maps  for  impulsively  started  drops  with  Red  =  1935.  The  density  ratios  are 
10  and  1.15  in  (a)  and  (b),  respectively. 

FIG.  14.  Eo  effect  on  drop  deformation:  pd/po  =  1-15,  Oho  =  0.05,  Ohd  =  0.0466.  The  drops 
are  accelerated  by  a  constant  body  force.  The  drop  evolutions  are  shown  for  four  Eo.  The  time  intervals 
between  successive  interfaces.  At*  are  (a)  2.5;  (b)  3.953;  (c)  5.590;  (d)  3.873.  The  grid  are  256  x  1280 
in  (b)  and  (c)  and  256  x  768  in  (a)  and  (d).  The  dashed  line  in  (a)  and  (d)  represents  the  actual  bottom 
boundary  of  the  computational  domain. 

FIG.  15.  Vorticity  contours  (left)  and  streamlines  (right)  for  the  drops  shown  in  Figure  14  (b)  and  (d). 
The  Eo  =  24  drop  is  plotted  at  t*  =  63.2  and  the  Eo  =  144  drop  is  plotted  at  t*  =  1 1.6. 

FIG.  16.  Non-dimensionalized  surface  area  plotted  versus  t*  for  the  drops  shown  in  Figure  14  (a)-(d). 
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FIG.  17.  Simulations  of  the  heat  transfer  during  drop  breakup  for  Eo  =  2.4, 24,  and  144.  Temperature 
field  and  drop  contour  are  plotted  for  three  drops  shown  in  Figure  14  (a),  (b),  and  (d)  at  selected  times. 
Oh  =  0.05,  Pr  =  1,  Pd/ Po  ~  1*15,  Pd/ Po  —  0.2,  Cp^f Cp^,  =  1,  and  k^/kg  —  0.2 

FIG.  18.  Average  temperature  plotted  versus  time  for  the  drops  shown  in  Figure  17. 
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